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ABSTRACT 



The cosmological many-body problem is effectively an infinite system of gravitation- 
ally interacting masses in an expanding universe. Despite the interactions' long-range 
(-^ nature, an analytical theory of statistical mechanics describes the spatial and veloc- 

ity distribution functions which arise in the quasi-equilibrium conditions that apply to 
O many cosmologies. Consequences of this theory agree well with the observed distribu- 

tion of galaxies. Further consequences such as thermodynamics provide insights into 
the physical properties of this system, including its robustness to mergers, and its tran- 
sition from a grand canonical ensemble to a collection of microcanonical ensembles with 
^ negative specific heat. 

1. Introduction 

(N 
O 

o> 

(One composed of many: Virgil, Moretum, 1, 104) 

> 

^ Imagine an expanding universe filled with objects moving under their mutual gravitational 

d forces. What decides the distribution of these objects in space, and their velocities, after the 

memory of their initial state has faded? Richard Bentley, England's leading 17th century classicist 
essentially posed the spatial part of this question to his friend Isaac Newton. Their ensuing letters 
may be found today in the library of Trinity College, Cambridge, and are discussed in some detail 



E pluribus unus 



elsewhere (Saslaw 2000), along with the subject's subsequent history. Briefly, Newton surmised 
that if the universe were finite the objects would eventually all fall together into one large cluster. 
But if the universe were infinite, they "could never convene into one mass; but some of it would 
convene into one mass and some into another, so as to make an infinite number of great masses, 
scattered at great distances from one to another throughout all that infinite space." And there the 
question rested for almost 300 years. 
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Today we restate this question by replacing "finite" and "infinite" by "static" and "expanding" 
and talking about galaxies instead of stars. The cosmological many-body problem appears to be 
a major component of the clustering of galaxies, although a detailed analysis of the astronomical 
observations also involves the roles of dark matter, galaxy formation and evolution, and perhaps 
dark energy. 



Early investigations of the gravitational many-body problem used a fluid approach when IJeans 



( 1902 1 explored the stability of a self-gravitating gas. In the cosmological problem, the fundamental 



particles are usually galaxies. Jeans' results described gravitational collapse, and gave a timescale, 
proportional to {Gp)~^^'^, for the collapse. 

Jeans' solutions were developed for a static universe, since expansion had not been discovered. 



However, Eddington (19301 showed that a static universe is unstable and a slight perturbation will 



cause it either to start expanding or collapsing. Bubble's ( 1929 1 discovery of expansion and the 



possible existence of dark energy (Riess et al. 1998; Perlmutter et al. 19991 provide a framework 



for the expanding universe and suggest that the universe will continue expanding so that structure 
formation and growth may eventually cease. This leads to interesting properties that we will explore 
below. 

In addition to expansion, the other major discovery that simplifies the cosmological many- 
body problem is that the two-point correlation function of galaxies, ^2{t), has a power law form 
^2(^") oc (Totsuji and Kihara 19691. This two-point galaxy correlation function is defined by 
the relation 

P{r\Ni = l)dr = 47rr^n{l + (,2{r))dr (1) 

for a statistically homogenous system of average number density n. P(r|A''i = 1) is the conditional 
probability given a galaxy in an infinitesimal volume at an arbitrary coordinate origin, that there 
is another galaxy at a distance r. Conventionally, spherical volumes are used although this can 
easily be generalized to volumes of arbitrary shape which provide further information. Two-point 
correlation functions therefore represent an excess over a random Poisson probability. In our case 
this excess is caused by the galaxies' mutual attraction. For positive 7, ^2('^) decreases significantly 
at large r. Observationally there is a correlation length, Ri, beyond which the correlation function 
decreases even faster than a power law, and may be neglected. Therefore over sufficiently large 
scales, the spatial distribution of galaxies is uncorrelated, and has a Poisson distribution, modified 
by its clustering on smaller scales. Modern surveys such as the 2DFGRS prawkins et al.||2003D have 
shown that 

/ \ -1.67 

" ( 5.05/. Mpc ) P> 

with Ri ~ 20h^^ Mpc and Hq = lOO/i where Hq is the Hubble constant. Amusingly, the exponent 
and amplitude in equation ^ are very close to those found in Totsuji and Kihara's analysis of the 
Lick survey, 40 years ago. 



This result implies that on sufficiently large scales, the average density and the gravitational 
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mean field is constant and isotropic. The resultant net force on a galaxy from the distant universe 
is negligible. On local scales, galaxies are subject to the long range effects of gravity from their 
neighbours within Ri. 



For an expanding universe, Saslaw and Fang (19961 showed that the expansion of the universe 



cancels the effect of the long range gravitational field for distances greater than Ri. This result 
is valid for Einstein-Friedmann models and for models that incorporate a cosmological constant, 
including ACDM cosmology. Because of expansion, the infinite range gravitational force has a finite 
effective range. Large-scale structures that form in the universe may be essentially stable as long 
as the universe is statistically homogenous at scales larger than Ri. 

These results simplify the problem. Since the expansion of the universe effectively cancels the 
long range gravitational field, we can calculate the potential on a galaxy by integrating over a finite 
region of space instead of over the entire universe. Although the problem still involves an infinite 
volume, the gravitational field now has a finite range and we can consider a finite sub-volume caused 
by this cancellation. 

Gravitational systems with more than two interacting particles are essentially unstable and 
not in equilibrium. While thermodynamics describes equilibrium systems, the cosmological case is 
characterised by quasi-equilibrium. This means that macroscopic quantities such as average tem- 
perature, pressure and density satisfy equilibrium relations whose variables change slowly compared 
with local relaxation timescales. For example, the average density of the universe(including dark 
matter) is approximately p = 3.5 x lO^'^mQMpc"^ ( Spergel et aL|[2007 |, and an average galaxy has 



a mass of approximately W^^niQ. Hence the dynamical timescale of the universe is ^"universe ~ 
Gyr, and there are about 0.35 galaxies per cubic megaparsec. The average peculiar velocity of a 
galaxy is ~ 1000 km s~^, or 1 Mpc Gyr^^. Therefore a cube with sides on the order of Ri will have 
about 3000 galaxies, and the time for a galaxy to cross the cube is about 20 Gyr, or 1.5 times the 
age of the universe. 

Local timescales are much shorter. A typical cluster similar to our local group of galaxies has 
a density(including dark matter) pic ~ 1-5 x lO^'^rnQMpc^^ . This is ~ 40 times greater than p, 
and hence the local group has a dynamical timescale tlg ~ ^O^^^^'^universe ~ ^ Gyr. Rich clusters 
which are much denser than the local group will have correspondingly shorter dynamical timescales. 
Most clusters have an average diameter of 2 ~ 4 Mpc, and hence the time taken to cross from one 
end to another is on the order of 3 Gyr. These numbers mean that "microscopic" perturbations on 
local scales will relax significantly faster than the "macroscopic" scale of a cube of side Ri, so that 
the system generally evolves from one equilibrium state to another in quasi-equilibrium. 



1.1. The GQED 



The preceding considerations frame the problem. We have a gravitational field with an effective 
range Ri in a system that is in quasi-equilibrium for which we need the counts-in-cells distribution 
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f{N), i.e. the probability that a given ceh located randomly in space contains galaxies. Thermo- 
dynamics and statistical mechanics provide two approaches to this problem. The thermodynamic 



solution was the first to be investigated (Saslaw and Hamilton 1984) and gives 



f{N-N,h) 



^^Ikl^ []V(l -h)+ Nh]""-' e-(^(i-^)+^^) 



(3) 



where N is the average number of galaxies in a cell, and 6 is a clustering parameter equal to the 
ratio between the correlation potential energy and twice the kinetic energy: 



W 
'2K 



boTiT 



1 + bonT 



(4) 



with < 6 < 1. Here T is the kinetic temperature of the system with kinetic energy K = 
3NT/2 (taking the Boltzmann constant = 1), and n = N /V is the average number of galaxies 
per unit volume where V is the cell volume. The distribution given by equation (|3]) describes the 
clustering of galaxies in quasi-equilibrium, and hence is known as the gravitational quasi-equilibrium 
distribution(GQED). 

The galaxy spatial distribution function f{N) is a simple but powerful statistic which char- 
acterises the locations of galaxies in space. It includes statistical information on voids and other 
underdense regions, on clusters of all shapes and sizes, on the probability of finding an arbitrary 
number of neighbours around randomly located positions, on counts of galaxies in cells of arbitrary 
shapes and sizes randomly located, and on galaxy correlation functions of all orders. These are just 



some of its representations (Saslaw 20001. Moreover it is also closely related to the distribution 



function of the peculiar velocities of galaxies around the Hubble flow (Saslaw et al. 1990 Leong 



and Saslaw 2004). 



While the form of b in the second equality of equation Q was originally taken as an ansatz by 



Saslaw and Hamilton (1984), a physical reason for the form was given by Saslaw and Fang (1996) 



using constraints from thermodynamics and the boundary conditions of the problem. Writing the 
correlation potential energy in terms of the two-point correlation function gives 



~ ~2K ~ 3r 
where m is the mass of an individual galaxy. 



V 



(5) 



The integral in equation ^ indicates that b depends on the shape and size of the cell, often 
taken to be spherical for simplicity. More generally, the shape of a cell affects the correlation 



potential energy so that (Saslaw 2000) 



Grin?ii? 
6NT 




V JV 



■dridT2. 



(6) 
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1.2. The Statistical Mechanical Approach 



The statistical mechanical approach to this problem was originated by Ahmad et al. ( 2002 1 who 
treated the cells as a grand canonical ensemble with the universe as an energy and particle reservoir 
and the galaxies as particles in the ensemble. The statistical mechanical results are essentially the 
same as those in the thermodynamic analysis of point particles, but here we use their generalization 
to extended particles. In this section we consider the simplest case with the following assumptions: 



• Galaxies behave like point masses and the two-galaxy interaction potential is 0(r) = —Gm? 



r . 



This will be generalised to include a softening parameter, e, at small scales in equation (23). 
But since the lower limits of the integrals in equations ([7|-([9]) converge uniformly in the limit 
€ = 0, we first use this limit to give the simple formulae in equations ( 10 )-(22 ) of the original 
thermodynamic results. 

• Only 2-body interactions are considered as the dominant type of interactions. 

• (j){r)T^^ is small. 

All three assumptions can be relaxed and we will discuss the resulting modifications in section [2| 

In order to compute the partition function, we note that since expansion cancels the long range 
mean-field force the integral over physical space has a finite cutoff and the integral 



Ri 



(7) 



exp [-0(r)r"^] = I exp d^r 
Jo I . 

for a galaxy mass potential has no singularities if the third assumption above holds. 

In order to compute the partition function, we start with the normalised phase space integral 

1 



Zn{T,V) 



exp 



1 /2vrmT\3^/2 




-1 



QNiT, V) 



(8) 



iV! V A2 J 

with A a normalisation constant and V the volume of the ensemble. From the observation that most 
gravitational interactions are dominated by two-body interactions, and by taking the first order 
expansion of the exponential, the potential energy term of the configuration integral Q]\f(T,V) 
becomes 



Qn{T,V) 



I exp[- (<^(ri,...,r;v))T-i] d'^r 



(ri,)T 



l<i<j<N 

n 

l<j<j<Af 



1 + 



Gm^ 



(9) 
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assuming that galaxies at large distances behave like point masses(although their finite extent is 
important in close interactions and mergers), that 2-body interactions dominate, and that (p{rij)T~^ 
is small. This gives the partition function for a canonical ensemble, the Helmholtz free energy, and 
the corresponding thermodynamic variables for entropy, pressure, internal energy and chemical 
potential: 

Z»(T,V)=^(^)"'V"(l-6)'-" (10) 
F = -TlnZK{T,V) = iVriii (^7^) + NTln (1 - b) - NT - ^-NTln (11) 



/dF\ NT 
U = F + TS = ^NT{l-2b) (14) 



T,V 

In these derivations the functional form of b occurs naturally as 

(3/2)G3m%r-3 
~ l + (3/2)G3m6nr-3' 



(16) 



Comparing the coefficients of equations (16) and (Q, we see that 

bo = {3/2)G^m^ (17) 

which quantitatively relates the correlation potential energy to the mass of an individual galaxy 
and con&ms the original ansatz in equation (Q. 

The distribution function is obtained by summing over all energy states for given N: 

_ e^^/^Z^(T,y) 



By using equations (10)-(15) and the relation between the grand canonical partition function Zq 
and the pressure equation of state, we then again obtain the distribution function given by equation 
(|3]). Figure [T] illustrates the counts in cells distribution for different values of b. For 6 = 0, there is 
no interaction and equation ([s]) becomes a Poisson distribution. 



- 7- 




1.3. The Peculiar Velocity Distribution 



The GQED also contains a consistent distribution function for the peculiar velocities of galax- 
ies. Since the partition function separates into kinetic and potential energy factors, the phase space 



distribution is separable such that (Saslaw et al. 1990) 



f{N,v)^ f{N)h{v). 



(19) 



To relate the counts- in-cells distribution to the peculiar velocity distribution, we need a relation 
between the number density N and the velocity v. This arises from the proportionality between 
the fluctuations in potential energy (due to correlations) over a given volume and the local kinetic 
energy fluctuations: 

GmN Q-^ = av'^ . (20) 

Here a is a constant of proportionality such that a = (1/r) (f^) ^ , v is the peculiar velocity of a 
galaxy, and r is the separation between two galaxies ( Leong and Saslawl|2004 ) . 

With the relation between and v in equation (20 1, the distribution in equation ^ can be 



transformed to a distribution in v (Saslaw et al. 1990): 



„, , , 20^/3(1-6) r ^. , , aiai-^-i _ 

1 \av^ + Ij 



(a/3(l-6)+af,.2)^^^ 



(21) 
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where (3 = (^2), ris the standard gamma function, and averages are taken over the grand canonical 
ensemble. Figure [2] illustrates the peculiar velocity distribution for different values of b. 

Since the proper motions of galaxies are too small to be observed, we can only measure the 
radial component of their peculiar velocity along our line of sight. We can generally write the 
velocity as a component parallel to our line of sight and a component perpendicular to our line 
of sight such that v'^ = v'^^ + v±- Then to obtain a form of the velocity distribution directly 
comparable with observations, we integrate over all perpendicular velocities to get the radial velocity 



distribution function( Inagaki et al. 1992) 



v± 



a/3(l -b) + ab (vf^ + vl^ 



r 



a[vf^+v^, 1+1 



— ab(v?i+v^, ) , 

-e V II ^ydv±. 



(22) 



Figures [9] and 



11 



below compare this with simulations and observations. 
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2. Modifications and the General Form of the GQED 

Various modifications and extensions exist for the GQED. The simplest case involves a mod- 
ification of the potential of a galaxy so that instead of treating galaxies as point masses, galaxies 



are treated as extended objects with a potential having the form ( Ahmad et al. 2002 1 

Grn^ 



(23) 



(r2+ £2)1/2 

where e is a softening parameter related to the radius of a galaxy. Such a modified potential is 
commonly used in A^-body simulations to model extended galaxies possibly with dark matter halos, 
and to avoid the singularity in the point mass potential when the separation between galaxies is 
small. The result of such a softening parameter is a modification h to such that 

^ {?,/2)G^m^C{e/Ri)nT-^ , . 

' l + (3/2)G3m6C(e/i?i)nr-3 ^ > 

where C, {e/Ri) is a term that depends on the interaction potential between a pair of galaxies. For 
X = e/Ri 

C {x) = \/l + x2 + In — — ^ . (25) 



In the case of a point mass, e ^ and C,{e/Ri) 1 so the partition function and its consequent 
thermodynamic properties converge to the solution for a point mass potential. Figure [3] illustrates 
C{e/Ri). 

To generalise this result, we can modify the point mass potential so that 

(j) = K{r) (26) 
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where K{r) is a dimensionless modification factor. It can then be shown that depends on the 
form of «;(r) such that — > 1 as ^(r) 1. This modification of the potential will only affect 6o- 
Therefore many other forms of the potential, such as multipole moments, are possible. 



2.1. Triplet Interaction 



Ahmad et al. (200661 investigated the contribution of the irreducible triplet interactions, which 



are the cases where three-body interactions occur very close together. In such a case, the modified 
partition function for a canonical ensemble of > 3 galaxies is 



1 



1 



N-l 



+ 



{N- l)(iV-2)4 



1 



(27) 



where we note that at large values of N, the {b/{l — b))^ term becomes very small compared to 
the (1/(1 — b))^~^ term. This indicates that the contribution from irreducible triplets is negligible 
compared to pairwise interactions for large N. 

From the partition function, the resulting distribution function taking into account irreducible 
triplets is 



f{N;N,b) 



N{1 - b) [N{1 -b) + Nb]^ V N^N^ ^L{N) 



-{N{l~bt)+Nbt) 



N\{1 + L{N)) 

where we write gn = 2{N — 1){N — 2)/9, and define L{N) and bt as follows: 



(28) 



L{N) 



aNb^{l - 6)^~^ for iV > 3 
for iV < 3 



(29) 



Af+3aivb^(l-b)^-3 
Af+AfaArb3(l-b)iv-3 



for iV > 3 
for < 3 



(30) 



Figure |4] compares the contribution of triplet interactions to the original GQED for different 
values of b. 



2.2. Two Species 



In order to help understand effects of a distribution of galaxy masses, Ahmad et al. (2006a 



derived the distribution function for a system consisting of two different galaxy masses and pointed 
out how this can be further generalised to a distribution of masses. We let the total number of 
galaxies in a cell be iV = A^^i + A''2 of which A'^i and N2 have individual masses mi and 1712 . 
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Fig. 4. — The GQED including irreducible triplets (equation 28) for N = 15(solid line). The dashed 
line, for comparison, is the GQED (equation [3| , without including irreducible triplets. 



By considering interactions between the different species, the partition function for a canonical 
ensemble of N galaxies becomes 



ZNiT, V) 



V 



N 



(27rmiT)3^i/2 (27rm2T)^^2/^ [l + bonT 



3N2/2 



-31^1-1 



[1 + bi2nT 



-31^2 



(31) 



where for the case of a Newtonian point mass potential, 60 = (3/2)(Gmf )^ and 612 = (3/2)(Gmim2)^. 
The distribution function is then 



f{N;N,b) 



m 

N{1 -b) + {m2/mifNb' 



N2 



1 — 6 + {m2/mi)^b 
where b is given as always by equations and (|5]) and B is given by 



^-(Nil~B)+NB) 



B 



{I + N2/N1) 



1 + 



{m2/miY{N2/Ni)' 
1 — 6 + {m2/mi)^b 



(32) 



(33) 



We compare the 2-species distribution function with the original GQED in figure [5} 
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Fig. 5. — The two species GQED from equation (32) for N = 15(solid line) with a number ratio of 
N2/N = 0.1 and mass ratio of 771-2/01,1 = 10. The dashed line, for comparison, is the GQED(equation 
3\ with mass mi and N = 15. 



2.3. Higher order expansions of exp(— (/>(r))/T 



A higher order expansion of the exponential in equation (|9|) is considered by Saslaw and Ahmad 
(2009 1. Again writing bo = (3/2)(Gm^)^, the canonical partition function is now 



Zn{T, V) 



Nl ( A2 



j [1 + bonT-\i + {bonT-^)\2 



(34) 



where and (2 are factors that arise from the second order expansion. In the case of point masses, 
Ci = 1 and C2 = 2/3. 



A quantity 6^ representing the modification of b by the second order expansion can be defined 



as 



6onr-3Ci + 2(6onr-3)2^2 



6(l-6)Ci+262C2 



1 + bonT~Hi + (6onr-3)2 (1 - b)^ + 6(1 - 6)Ci + b^C2 
The distribution function in this case is 



(35) 



f{N;N,b) 



N{1 - b) 

m 



N{l-b)+NbQi + ^X2 



N{l-b)' 

(l-6) + feCi+(S)C^ 



N-l 



^-{N{l-b^)+NK) 
(l-6)+Kl+ (l5)C2' 



(36) 
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Normally the effects of the high order terms in this expansion are small, thereby confirming 
the essential form of f{N) in equation Q and its consequences. 



3. Properties of the GQED 



The GQED has a number of properties that are discussed in detail in Saslaw (20001. We 
summarise some of them here. 



3.1. Correlation Functions 



For the simplest case, the two-point correlation function enters f{N) through (joj) since b 
depends on the volume integral of ^2- However, in general the form of the GQED also contains 
information about the volume integrals of all the correlation functions where the average volume 
integral is defined as 



N 



1 



V 



^Ar(ri, . . . ,rAr)(i^ri ...tn- 



(37) 



For the general A^-point correlation function, Zhan and Dyer (19891 derived a relation between 
A'^ and b. The general form for is 



N 



1-b 

bN^ 



M=N 



M 



M-l 



{M -N)\ 



For = 1, 2 and 3 this gives 



^3 



b 2-6 
= (1-6)2 ]v 

62 9 - 86 + 262 



(1-6)^ 



(38) 

(39) 
(40) 

(41) 



Average moments of other quantities, including the thermodynamic variables, have also been de- 
rived. 



3.2. Specific Heat 



From the internal energy given by equation ( 14 ), using equation (Q the specific heat per galaxy 
at constant volume is 

=^(1 +46 - 662) (42) 



Cv 



N dT 
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Fig. 6. — The specific heat per galaxy at constant volume for varying b from equation (42) 



which at 6 = 0, is 3/2 representing the monatomic ideal gas. As galaxies cluster, binaries start to 
dominate and Cy reaches its maximuiTi of 5/2 at b — 1/3. At a. critical value of bcrit 

such that 



crit 



2 + VIO 
6 



0.8604 



(43) 



Cy = 0. As 6 increases beyond bcrit-, Cy decreases until Cy reaches —3/2 at 6 = 1 representing 
the specific heat of a fully virialized cluster. The transition from positive to negative specific heat, 
illustrated in figure [6] occurs at a rather high value of 6, and does not involve a discontinuity, unlike 
some laboratory systems. 

While negative specific heat is frequently encountered in the gravitational context, the system 
we are discussing is nominally a grand canonical ensemble. Grand canonical ensembles can only 
have positive specific heat( Tblman|1938 equation 141.12). However for a microcanonical ensemble, 
the specific heat can be negative! Thirring||1970p (For a review see Dauxois et al. (2002 1). 



To resolve this "paradox", we consider clustering for increasing b. As b increases, galaxies 
become increasingly clustered and large voids are created. These voids have the effect of insulating 
clusters from the energy and galaxy reservoir that is the rest of the universe. The net effect is 
that the universe breaks up into a collection of galaxy clusters. Each cluster is a microcanonical 
ensemble that does not exchange significant energy or galaxies with the rest of the universe. 

The measured value of b for the local universe of redshifts z < 0.1 is about 0.86 for large 8° 
angular cells ( SivakofF and Saslaw||200"5 1 which indicates that b ~ bcrit and suggests that clustering 
is fairly advanced in the local universe. The average matter density(including dark matter) of the 



universe is 3.5 x lO^*^m0Mpc ^ (Spergel et al. 2007 1 . For comparison, the mass and radius of a 
typical rich cluster are about lO^^m© and 2 Mpc, making it approximately 1000 times denser than 
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the average universe. To find a lower bound to the inter-cluster spacing, consider a simple model in 
which rich clusters are each surrounded by a void where no galaxies exist. Then each void would be 
approximately a sphere that would contain roughly lO^^m© if the matter were initially distributed 
uniformly throughout the universe. The radius of such a sphere is approximately 19 Mpc, which 
suggests that the average massive cluster separation is at least 19 Mpc. 

For an average galaxy with a peculiar velocity of about 1000 km/s to move from one cluster 
to another over a distance of 19 Mpc would require about 19 Gyr, which is longer than the age 
of the universe. The long transit times involved suggest that such clusters are no longer likely 
to exchange galaxies. At separations of 19 Mpc, the interaction potential energy of a pair of 
galaxies of ~ IO^^ttiq each is approximately 4 orders of magnitude lower than the kinetic energy 
of a single galaxy with a peculiar velocity of 1000 km/s which indicates that at large values of b, 
clusters would exchange energy more effectively by exchanging galaxies rather than by long range 
gravitational interactions. But galaxy exchange now takes too long. Hence at 6 > bcrit, clusters are 



approximately microcanonical ensembles that are likely to be virialized(Leong and Saslaw 2004 1 
and have a negative heat capacity. When most galaxies are bound to a cluster in the case of large 
b, the total heat capacity is negative, and hence the specific heat per galaxy is negative. 

To see this, we can use the multiplicity function derived from equation ([s]) which gives the 



probability that a physical cluster contains galaxies (see Saslaw 2000 equation 28.97): 



ri{N, b) = ^^^^ ' e-^^ for A = 1, 2, 3, . . . (44) 



which is a truncated Borel distribution. From equation(42 ), the total heat capacity of these clusters 
is therefore 

o oo _ ^ 

- (1 + 46 - 662) ^ NNor]{N, 6) = - AA^o (l + 46 - 66^) (45) 

N=l 

where Aq is the total number of clusters having an average number A of galaxies per cluster. The 
multiplicity function sums to unity and we can divide the total heat capacity by A^Aq to obtain the 
average specific heat of the collection of clusters, each of which is a microcanonical ensemble. This 



gives the same result as equation ( 42 1 , so the effective specific heat of this ensemble also becomes 



negative for b > bcrit in equation (43). 



When 6^1, /(A) for A > 0. The void probability goes to 1 and all galaxies will be 
bound to a single cluster that can be represented as an isothermal sphere( Baumann et al.|2003 |. In 



that limit, the energy and galaxy reservoir is effectively depleted, and the entire system effectively 
becomes a single microcanonical ensemble. 

This transition, as 6 increases, from a single grand canonical ensemble where the average 
cluster has a positive specific heat per galaxy through a collection of microcanonical ensembles 
with negative specific heat per galaxy to a single virialized microcanonical ensemble is an effect 



of increasing gravitational inhomogeneity. It is rather remarkable that equation (42) can describe 
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this entire ensemble transition in a smooth manner. Perhaps the reason is that properties of the 
system for b > bent are essentially determined by its properties for b < bcrit- 



3.3. Robustness to Mergers 

Galaxies are known to interact and merge with each other. Over time, we can expect the 
average number of galaxies in a given comoving volume to decrease due to mergers, and hence 
galaxy mergers will change the distribution of galaxies. We can classify galaxy mergers by the 
masses of the progenitors mi and m2. Minor mergers are mergers where one galaxy is much less 
massive than another galaxy such that mi <^ m2 or m2 <^ mi. In such mergers, the resulting 
galaxy will have a position that is close to the centre of mass of the larger of the two galaxies 
and the only change to the distribution is the reduction of the number of low mass galaxies. The 
other class of mergers are major mergers where both progenitors are of comparable mass such that 
mi ~ m2, and the position of the resulting galaxy is approximately at the midpoint between the 
two. Major mergers will change the spatial distribution of galaxies more drastically. The average 
number density N decreases with time and the average mass of each galaxy increases with time. 
This affects the rate at which b changes. 

To follow the change in the distribution of galaxies due to mergers, we consider the distribution 
of nearest neighbour galaxy pairs whose members are separated by a distance 2r. We denote the 
separation between the midpoints of two such pairs by R. Hence the interaction between two 
galaxy pairs, each with separations 2ri and 2r2, is given by 

Gm^ ( |R| IRI IRI |R| 

+ + ,„ ' . r + 



|R| V|2ri| |2r2| |R + ri+r2| |R - ri + raj 

R R \ / . V 

+ |R + r.-r.| + |R-r.-rJ 

where the term in brackets can be viewed as a modification to the point mass potential. By averaging 



over possible values of ri and r2, we can use equation (46) to define a modification k(|R|, (r)) to 
the potential where |R| is the separation between two pairs. From section |2] The modification due 
to such a potential will enter into the form of 6o and hence the resulting distribution will also follow 
the GQED. By assuming that all galaxies will eventually merge, we note that when all galaxies 
have merged once, the spatial distribution of the resulting galaxies also follows the GQED. Hence 
the distribution of galaxies is robust to mergers. 



3.4. The Time Evolution of b 



Since there is no final equilibrium state for a classical infinite gravitating system, its clustering 
increases as described by the increase of b with time. From the thermodynamic variables in section 
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[T] we can write (cf. [Saslaw (1992)) 



b = bo{N)PT-*. (47) 
Then by taking the expansion of the universe to be an adiabatic process, 

= dU + PdV = -{1 - 2b) TdN\T,P + NdT\j^^p - SNTdb + NT{1 - b) — . (48) 



If galaxy mergers do not occur, dN = 0. With the relation dV/V = 3da/a, we arrive at a 
relation relating the scale length, o, of the universe with b and N such that 



l + 6bdb 1-6 



86 da 

Integrating equation (49) we get( Saslawl|1992 ) 

51/8 



(1 - 6)7/8 a, 

where a^: is a constant of integration given by the initial state. 

Galaxies however will merge, and hence we can expect that dN 7^ 0. This gives 



(49) 



(50) 



1 + 66 d6 1-6 l-2b dN 
+ 



86 da 



2N da 



(51) 



From the definition of 6 given in equation (|4j), we note that 6 depends on n and hence and V. 
Denoting rn as the average mass of a galaxy, by conservation of mass we have fnN = constant for 
an ensemble of comoving volumes. This gives 



db 

T 



^ N 



(52) 



Then from equations (51 ) and (52), the rate of change of 6 is given by 

1 + 66 1-26 



86 



+ 



126C4e/i?i) J da a 



db 



1 



where is a term that is given by 



Ue/Ri 



1 + 



1 6 C'{^/Ri 
18 Ri Q{e/Ri 



(53) 



(54) 



for the case of the isothermal halo with Q{e/Ri) given by equation (25). Here C,'{e/Ri) is the 



derivative of Q{e/Ri) with respect to e/Ri- In this case is a factor of order unity that varies 
between 1 and 17/18. 
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In the limit of a small halo, e/Ri 
become 



and C^(e/i?i) 
55/24 



a 
a* 



1 and equation (53) integrates to 



(55) 



(1 - 6)19/24 

In both the non-merging and merging cases, we note that b increases as a increases with b increasing 
faster in the case with no mergers than in the case with mergers. We compare the two cases in 
figure [7} 



4. Simulations and Observations 

The results of the GQED have been compared with A^-body simulations and various observa- 
tions, and there is very good agreement between predictions, simulations and observations without 
having to introduce additional(e.g. non-gravitational) parameters. 



4.1. A/^-body Simulations 



A series of computer simulations (Itoh et al. 1988 1990 Inagaki et al. 1992 Itoh et al. 19931 
examined the various aspects of the GQED in a comoving volume with N between 4000 and 10000. 
The spatial distribution and peculiar velocity distributions were examined for cases where galaxies 
had the same mass and for various mass spectra. The primary results of these simulations were 
that the distributions of galaxies follow the GQED very well for identical masses, and even better 
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Fig. 8. — fv{N) for a single simulation of JIq = 1 and cold initial conditions for the expansion 



factor a/aQ = 7.86 from Itoh et al. (1988). Various cell radii R between 0.1 to 0.4 are plotted. The 



solid lines are from the simulation and the dashed lines are from equation (|3]). 



for a more realistic case where galaxies have different masses. Figure [8] shows a typical example for 
identical masses and figure |9] shows examples of velocity distributions. 



4.2. Observations of fv{N) 

There have been many observations of fv{N) using various galaxy catalogs. Some recent 
ones (which also give earlier references) include 



1. An analysis of 2-dimensional angular cells (Sivakoff and Saslaw 20051, figure 10 using data 



from the 2MASS survey, an all-sky survey in the infrared with up to 439754 galaxies. A good 
fit to predictions was found, with, b approa/Cliing h^j-n at large cell sizes. 



2. A 3-dimensional analysis of the Pisces- Perseus super cluster with 4501 galaxies (Saslaw and 
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Fig. 9. — The radial velocity distribution f{vr) for a single simulation of = 1.0 from Inagaki et al. 



(19921. The histogram is data from the simulation. Various expansion factors a/ao are plotted. 
The solid lines are from equation (22) and the dashed lines are Maxwell-Boltzmann distributions 
with the same (v^). 
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Fig. 10. — The angular projected GQED fe{N) for square cells with different values of 9 from 



Sivakoff and Saslaw (2005). The solid histogram is data from observations. The dashed curve is 
from equation (|3|. Values of h and N were found directly from the observations. 



Haque-CopiTah||l998[ ) found N = 3.38 and h = 0.80 for cells of lO/i"^ Mpc where h is the 
reduced Hubble constant such that h = //q/IOO. The value was found to be 0.089, 
indicating a relatively good fit to the GQED. 



3. At high redshifts, Rahmani et al. (2009) found that the projected spatial distribution follows 



the GQED out to a redshift of z = 1.5 using the GOODS survey catalog. Due to the small 
sample size and limited sky coverage, large differences between the North and South fields in 
the GOODS catalog are evident. Nevertheless, there is a good agreement between the form 
of the GQED and the observed spatial distribution of galaxies even at these high redshifts. 



4.3. Observations of f{v) 



The radial velocity distribution was analysed observationally ( Raychaudhury and Saslaw|1996 ) 
using the Matthewson catalog of spiral galaxies. Since radial velocities of galaxies require a sec- 
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Fig. 11. — The observed radial velocity distribution f{vr) from Raychaudhury and Saslaw (19961 



The histogram is data from observations. The solid lines are from equation (22) and the dotted 



lines are best fit Maxwell-Boltzmann distributions. The left panel shows the radial velocity dis- 
tribution corrected for a bulk motion of 599 km/s while the right panel shows the radial velocity 
distribution corrected for an extra Hubble expansion with an effective local Hubble constant of 
Ho = 92 km/s/Mpc. 



ondary distance indicator and are not easily obtained, the catalog is relatively small with 1353 
galaxies. Despite the small sample, a relatively good fit to the radial velocity distribution of equa- 



tion (22) was found as illustrated in figure 11 



5. Conclusion 

Although the cosmological many-body problem is essentially a long range gravitating system, 
the expansion of the universe cancels the long-range gravitational mean field so that interactions 
effectively have a finite range. Using this property and its consequences, we have described the 
spatial and velocity distributions of particles in such a system. Its fundamental particles are 
galaxies. In the simplest case the spatial distribution of galaxies can be characterised by just the 
average number of galaxies in a cell, N, and a clustering parameter b, which is essentially the ratio 
of interaction and kinetic energies. This assumes that galaxies are equal point masses, and that 
pairwise interactions are the dominant form of interactions. Relaxing these assumptions introduces 
minor higher order corrections into b. However the form of the spatial distribution does not change, 
suggesting it is very robust to a wide range of physical conditions. 

The velocity distribution function can be derived consistently from the spatial distribution 
which can also be used to obtain the volume integrals of the A'^-point correlation functions. From 
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the thermodynamic quantities the time evolution of b can be derived, describing the rate at which 
clusters form. Towards increasing b, the grand canonical ensemble gradually breaks up into a 
collection of microcanonical ensembles, each of which are virialized clusters that individually have 
negative specific heat such that the average specific heat goes negative when b is greater than a 
critical value bcrit = 0.8604. 

Finally, comparisons with A'^-body simulations and observations show that the GQED does 
indeed describe the physical world, at both low and high redshifts and for both spatial and velocity 
distributions. 

Acknowledgments: We are both grateful to the organisers of the Les Houches Summer 
School for an environment of stimulating discussion. We particularly thank Thierry Dauxois and 
Michael Kastner for their comments on negative specific heat. 

REFERENCES 

Ahmad F., Malik M. A., and Masood S. (2006a). International Journal of Modem Physics D, 15, 
1267-1282. 

Ahmad F., Saslaw W. C, and Bhat N. I. (2002). ApJ, 571, 576-584. 
Ahmad F., Saslaw W. C, and Malik M. A. (20066). ApJ, 645, 940-949. 
Baumann D., Leong B., and Saslaw W. C. (2003). MNRAS, 345, 552-560. 

Dauxois T., Ruffo S., Arimondo E., and Wilkens M. (2002). In Dynamics and Thermodynam- 
ics of Systems with Long-Range Interactions (ed. T. Dauxois, S. Ruffo, E. Arimondo, and 
M. Wilkens), Volume 602, Lecture Notes in Physics, Berlin Springer Verlag, pp. 1-19. 

Eddington A. S. (1930). MNRAS, 90, 668-678. 

Hawkins E. et al. (2003). MNRAS, 346, 78-96. 

Hubble E. (1929). Proceedings of the National Academy of Science, 15, 168-173. 
Inagaki S., Itoh M., and Saslaw W. C. (1992). ApJ, 386, 9-18. 
Itoh M., Inagaki S., and Saslaw W. C. (1988). ApJ, 331, 45-63. 
Itoh M., Inagaki S., and Saslaw W. C. (1990). ApJ, 356, 315-331. 
Itoh M., Inagaki S., and Saslaw W. C. (1993). ApJ, 403, 476-496. 

Jeans J. H. (1902). Royal Society of London Philosophical Transactions Series A, 199, 1-53. 
Leong B. and Saslaw W. C. (2004). ApJ, 608, 636-646. 



- 24 - 



Perlmutter S. et al. (1999). ApJ , 517, 565-586. 

Rahmani H., Saslaw W. C, and Tavasoli S. (2009). ApJ . In press. 

Raychaudhury S. and Saslaw W. C. (1996). ApJ, 461, 514-524. 

Riess A. G. et al. (1998). AJ, 116, 1009-1038. 

Saslaw, W. C. (1992). ApJ, 391, 423-428. 

Saslaw W. C. (2000). The distribution of the galaxies : gravitational clustering in cosmology. 
Cambridge, U.K. ; New York : Cambridge University Press, 2000. 

Saslaw W. C. and Ahmad F. (2009). in preparation. 

Saslaw W. C, Chitre S. M., Itoh M., and Inagaki S. (1990). ApJ, 365, 419-431. 

Saslaw W. C. and Fang F. (1996). ApJ, 460, 16-27. 

Saslaw W. C. and Hamilton A. J. S. (1984). ApJ, 276, 13-25. 

Saslaw W. C. and Haque-Copilah S. (1998). ApJ, 509, 595-607. 

SivakofF G. R. and Saslaw W. C. (2005). ApJ, 626, 795-808. 

Spergel D. N. et al. (2007). ApJS, 170, 377-408. 

Thirring W. (1970). Zeitschrift fur Physik, 235, 339-352. 

Tolman R. C. (1938). The Principles of Statistical Mechanics. Oxford, U.K. : The Clarendon 
Press. 

Totsuji H. and Kihara T. (1969). PAS J, 21, 221-229. 
Zhan Y. and Dyer C. C. (1989). ApJ, 343, 107-112. 



This preprint was prepared with the AAS IM^jX macros v5.2. 



